subroutine solve_u_p_diedai(u1,p1,rho1,u2,p2,rho2,u,p)
implicit none
   
    real(8) u1,p1,rho1,u2,p2,rho2  !init
    real(8) u,p,old_p  !solver
    real(8) Fp,f1,f2
    real(8) ebuso
    real(8) a1,a2,c1,c2   !c is sound velocity
    real(8),parameter:: gamma=1.4d0
    real(8) b_v!between var
    
    ebuso=1e-10
    c1=dsqrt(gamma*p1/rho1)
    c2=dsqrt(gamma*p2/rho2)
    a1=rho1*c1;a2=rho2*c2
    old_p=0
    do 
        p=a1*p2+a2*p1-a1*a2*(u2-u1)
        p=p/(a1+a2)
        p=max(p,ebuso)
        if(dabs(p-old_p)<ebuso)then
            exit
        endif
        old_p=p
        call f_function(f1,p,p1,rho1)
        call f_function(f2,p,p2,rho2)
        Fp=f1+f2
        if(p<=ebuso)then
            exit
        endif
        if(p<p1)then
            b_v=(gamma-1.d0)/(2.d0*gamma)
            a1=rho1*c1*(gamma-1.d0)*(1.d0-p/p1)
            a1=a1/(1.d0-(p/p1)**b_v)
            a1=a1/(2.d0*gamma)
        else
            a1=(gamma+1.d0)/(2.d0*gamma)
            a1=a1*(p/p1-1.d0)+1.d0
            a1=dsqrt(a1)
            a1=rho1*c1*a1                
        endif
        if(p<p2)then
            b_v=(gamma-1.d0)/(2.d0*gamma)
            a2=rho2*c2*(gamma-1.d0)*(1.d0-p/p2)
            a2=a2/(1.d0-(p/p2)**b_v)
            a2=a2/(2.d0*gamma)
        else
            a2=(gamma+1.d0)/(2.d0*gamma)
            a2=a2*(p/p2-1.d0)+1.d0
            a2=dsqrt(a2)
            a2=rho2*c2*a2           
        endif
    enddo
    u=(a1*u1+a2*u2+p1-p2)/(a1+a2)

    
    
endsubroutine